В данном задании вам необходимо реализовать алгоритм кластеризации Partition Around Medoids.
Два возможных варианта реализации:
Нужно написать функцию, которая принимает на вход несколько параметров и возвращает также несколько значений.
Параметры функции:
pdist());Возвращаемые значения:
По аналогии с классами в scikit-learn, нужно реализовать класс, наследуемый от Base Estimator.
Подробнее про реализацию своих моделей в scikit-learn: here.
Параметры:
pdist());Методы:
Атрибуты:
Note 1: Параметры max_iter и tol должны иметь дефолтные значения.
Note 2: Функции для вычисления расстояний самим реализовывать не нужно, используйте pdist().
Также необходимо написать документацию к функции/методу: описать формат входных данных (параметров) и возвращаемых значений, особенности работы функции и детали реализации алгоритма. В качестве образца можно взять документацию к функциям/методам, которые рассматривали на занятиях.
Наивная реализация алгоритма будет работать довольно медленно - это нормально. Будет плюсом (но не является обязательным), если вы попытаетесь оптимизировать ваш код. Можете указать все ваши решения для оптимизации в документации.
import random
import pandas as pd
import pylab as pl
import numpy as np
import scipy.spatial as ss
import sklearn.cluster as sc
import sklearn.manifold as sm
import sklearn.datasets as ds
import sklearn.metrics as smt
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, cdist
from sklearn.base import BaseEstimator
from itertools import cycle
from math import hypot
def medoid_distribution(med_ind, data):
dist_vector = np.zeros(len(data)) # ближайшая медоида для каждого объекта
dist_matr = np.zeros((len(data), len(med_ind))) # матрица расстояний
target = np.Inf # значение целевой функции
min_dists = np.zeros(len(data))
for idx, p in enumerate(data):
''' поиск для каждого объекта данных ближайшего к нему медоида,
приписание объекта к этому медоиду: "в матрице расстояний" индексу данного объекта
присваивается индекс медоида'''
dist_matr[idx] = [np.linalg.norm(p - data[q]) for q in med_ind]
dist_vector[idx] = med_ind[np.asarray(dist_matr[idx]).argmin()]
min_dists[idx] = min(dist_matr[idx])
target = sum(min_dists)
return dist_vector, dist_matr, min_dists, target
def medoid_remake(med_ind, dist_vector, data):
new_meds = []
for m in med_ind:
clust = np.where(dist_vector == m) # индекс тех эл-в в данных, медоид которых равен m
clust = clust[0] # функция np.where() возвращает tuple из 1 эл-та, содержащего array
distances = [] # здесь будут суммы р-й каждой точки до всех остальных в этом кластере
for i in clust:
distances.append(sum([np.linalg.norm(data[i] - data[q]) for q in clust if q != i]))
new_meds.append(clust[np.asarray(distances).argmin()])
return new_meds
class PAM(BaseEstimator):
def __init__(self, k=3, metric="euclidean", max_iter = 300, tol=0.001):
self.k = k
self.metric = metric
self.max_iter = max_iter
self.tol = tol
def fit(self, dt):
# преобразование данных из dataframe, если необходимо в массив списков характеристик каждого объекта
if type(dt) == pd.core.frame.DataFrame:
data = dt[col].values.tolist()
else:
data = dt
# Шаг 1
med_ind = set() # рандомная генерация индексов медоидов
while len(med_ind) != self.k:
med_ind.add(random.randint(0, len(data) - 1))
med_ind = list(med_ind)
# Шаг 2 матрица расстояний
dist_vector, dist_matr, min_dists, TARGET_INIT = medoid_distribution(med_ind, data)
# TARGET_INIT - начальное значение целевой ф-и, с которым будем сравнивать её изменения
# Шаг 3 пересчёт расстояний и медоид
for i in range(self.max_iter):
med_ind = medoid_remake(med_ind, dist_vector, data)
dist_vector, dist_matr, min_dists, target = medoid_distribution(med_ind, data)
if TARGET_INIT - target < self.tol:
break
self.inertia_ = target
self.medoids_ = med_ind
self.labels_ = dist_vector
#return self.inertia_, self.medoids_, self.labels_
# pam = PAM(2)
# pam.fit(data_pca[: 50])
# pam.labels_
# pam.medoids_
# pam.inertia_
В рамках данной лабораторной работы вам предлагается проанализировать набор данных по различным городам США. Каждый город характеризуется следующими признаками:
import pandas as pd
pd.set_option('display.max_rows', 20)
pd.set_option('display.max_colwidth', 1000)
data_desc = pd.read_csv('Data_Description.txt', sep=':')
data_desc
Housing и Crime - наоборот.Population- статистический признак, не имеющий интерпретации как “лучше-хуже”.Place - уникальный идентификатор объекта (города), он не должен использоваться при кластеризации.Longitude и Latitude. Их также не следует использовать при кластеризации данных.df = pd.read_csv('Data.txt', sep=' ').reset_index()
df
0. Выполните необходимую предобработку данных. Перед кластеризацией исключите из данных признаки Place, Long и Lat.
data = df.drop(["Long", "Lat", "Place", "Pop"], axis=1)
data['HousingCost'] = -1 * data['HousingCost'] # *(-1), тогда значения будут интерпретироваться как для всех остальных
data['Crime'] = -1 * data['Crime']
data
from sklearn.preprocessing import StandardScaler
data_scaled = StandardScaler().fit_transform(data)
from sklearn.cluster import *
import scipy.cluster.hierarchy as sch
import mpl_toolkits.mplot3d.axes3d as p3
import plotly.express as px
from sklearn.decomposition import PCA
pca = PCA(n_components=0.9, random_state=321)
data_pca = pca.fit_transform(data_scaled)
from itertools import combinations
col = list(data.columns)
dendrogram = sch.dendrogram(sch.linkage(data_pca, method='ward'))
dendrogram = sch.dendrogram(sch.linkage(data_pca, method='complete'))
for i in range(10, 34):
ward = AgglomerativeClustering(n_clusters=None, linkage='ward', distance_threshold=i).fit(data_pca)
print(*set(ward.labels_))
import random as rnd
# Визуализация для Climate и HousingCost (первые два столбца) разделения на 2-9 кластеров, linkage = 'ward'
for i in range(2, 10):
ward = AgglomerativeClustering(n_clusters=i, linkage='ward').fit(data_pca)
label = ward.labels_
fig = plt.figure(figsize=(8,8),)
for k in range(i):
x = rnd.randint(2, 255)
y = rnd.randint(2, 255)
z= rnd.randint(2, 255)
plt.scatter(data_pca[label==k, 0], data_pca[label==k, 1], s=50, marker='o',
color='#{:02x}{:02x}{:02x}'.format(x,y,z))
plt.show()
# Визуализация для Climate и HousingCost (первые два столбца) разделения на 2-9 кластеров, linkage = 'complete'
for i in range(2, 10):
ward = AgglomerativeClustering(n_clusters=i, linkage='complete').fit(data_pca)
label = ward.labels_
fig = plt.figure(figsize=(8,8),)
for k in range(i):
x = rnd.randint(2, 255)
y = rnd.randint(2, 255)
z= rnd.randint(2, 255)
plt.scatter(data_pca[label==k, 0], data_pca[label==k, 1], s=50, marker='o',
color='#{:02x}{:02x}{:02x}'.format(x,y,z))
plt.show()
for i in range(2, 10):
label = AgglomerativeClustering(n_clusters=i, linkage='ward').fit(data_pca).labels_
fig = plt.figure()
ax = p3.Axes3D(fig)
ax.view_init(7, -80)
for l in np.unique(label):
ax.scatter(data_pca[label == l, 0], data_pca[label == l, 1], data_pca[label == l, 2],
color=plt.cm.jet(float(l) / np.max(label + 1)),
s=20, edgecolor='k')
for k in range(2, 6):
fig = plt.figure(figsize=(30,10))
label = AgglomerativeClustering(n_clusters=k, linkage='ward').fit(data_pca).labels_
j = 0
print(f"{k} clusters")
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
plt.show()
from sklearn.neighbors import NearestNeighbors
for i in range(3, 10):
nbrs = NearestNeighbors(n_neighbors=i, algorithm='ball_tree').fit(data_pca)
distances, indices = nbrs.kneighbors(data_pca)
plt.plot(sorted(distances[:, -1]))
plt.ylabel(f"k-distances for {i} neighbours")
plt.grid(True)
plt.show()
# для каждой точки высчитывается среднее расстояние до n ближайших соседей
def k_neib(p, data, n):
return sum(sorted([np.linalg.norm(p - q) for q in data])[: n]) / n
# для каждой точки в датасете высчитывается среднее р-е до n ближайших соседей
def all_knn(min_pts):
all_knn = []
for i in data_pca:
all_knn.append(k_neib(i, data_pca, min_pts))
return sorted(all_knn)
# функция all_knn(min_pts) возвращает отсортированный массив средних расстояний каждой точки до её n ближайших соседей
for i in range(3, 10):
nbrs = all_knn(i)
plt.plot(nbrs)
plt.ylabel(f"k-distances for {i} neighbours")
plt.grid(True)
plt.show()
db = DBSCAN(eps=0.9, min_samples=4).fit(data_pca)
label = db.labels_
#label
from itertools import combinations
col = list(data.columns)
# визуализация DBSCAN для каждой пары признаков при eps=0.9 - 6 кластеров - слишком много
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
db = DBSCAN(eps=1.25, min_samples=4).fit(data_pca)
label = db.labels_
#label
# визуализация DBSCAN для каждой пары признаков при eps=1.25 - тоже много
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
db = DBSCAN(eps=2, min_samples=4).fit(data_pca)
label = db.labels_
#label
# визуализация DBSCAN для каждой пары признаков при eps=2
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
db = DBSCAN(eps=4, min_samples=4).fit(data_pca)
label = db.labels_
#label
# визуализация DBSCAN для каждой пары признаков при eps=3 - не все выбросы считаются таковыми
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
# трёхмерное изображение по тройкам признаков
j = 0
for i in combinations(range(7), 3):
label = DBSCAN(eps=0.9, min_samples=4).fit(data_pca).labels_
fig = plt.figure()
ax = p3.Axes3D(fig)
ax.view_init(7, -80)
for l in np.unique(label):
ax.scatter(data_pca[label == l, i[0]], data_pca[label == l, i[1]], data_pca[label == l, i[2]],
color=plt.cm.jet(float(l) / np.max(label + 1)),
s=20, edgecolor='k')
plt.legend([data_desc['Attribute'][i[0]], data_desc['Attribute'][i[1]], data_desc['Attribute'][i[2]]])
j += 1
# подбор параметров методом локтя и силуэта
from sklearn.metrics import silhouette_score
km_scores= []
km_silhouette = []
vmeasure_score =[]
db_score = []
for i in range(2,10):
km = KMeans(n_clusters=i, random_state=0).fit(data_pca)
preds = km.predict(data_pca)
km_scores.append(-km.score(data_pca))
silhouette = silhouette_score(data_pca, preds)
km_silhouette.append(silhouette)
plt.figure(figsize=(30, 4))
plt.subplot(1, 2, 1)
plt.title("Elbow method\n", fontsize=15)
plt.scatter(x=[i for i in range(2, 10)], y=km_scores, s=150, edgecolor='k')
plt.grid(True)
plt.xlabel("Number of clusters", fontsize=15)
plt.ylabel("K-means score", fontsize=15)
plt.xticks([i for i in range(2, 10)], fontsize=15)
plt.yticks(fontsize=15)
plt.subplot(1, 2, 2)
plt.title("Silhouette coefficient method \n", fontsize=15)
plt.scatter(x=[i for i in range(2,10)], y=km_silhouette, s=150, edgecolor='k')
plt.grid(True)
plt.xlabel("Number of clusters", fontsize=15)
plt.ylabel("Silhouette score", fontsize=15)
plt.xticks([i for i in range(2,10)], fontsize=15)
plt.yticks(fontsize=15)
plt.show()
km = KMeans(n_clusters=2, random_state=123).fit(data_pca)
label = km.labels_
#label
km.cluster_centers_
for k in range(2, 6):
fig = plt.figure(figsize=(30,10))
label = KMeans(n_clusters=k, random_state=123).fit(data_pca).labels_
j = 0
print(f"{k} clusters")
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
plt.show()
from datetime import datetime
pam = PAM(2)
start_time = datetime.now()
pam.fit(data_pca)
delta = datetime.now() - start_time
pam.labels_
print(f"Self-written PAM algorithm work time: {delta.seconds} sec {delta.microseconds} microsec")
# трёхмерное изображение по тройкам признаков для n=2
j = 0
for i in combinations(range(7), 3):
fig = plt.figure()
ax = p3.Axes3D(fig)
ax.view_init(7, -80)
for l in np.unique(label):
ax.scatter(data_pca[label == l, i[0]], data_pca[label == l, i[1]], data_pca[label == l, i[2]],
color=plt.cm.jet(float(l) / np.max(label + 1)),
s=20, edgecolor='k')
plt.legend([data_desc['Attribute'][i[0]], data_desc['Attribute'][i[1]], data_desc['Attribute'][i[2]]])
j += 1
for k in range(2, 6):
fig = plt.figure(figsize=(30,10))
pam = PAM(k)
pam.fit(data_pca)
label = pam.labels_
j = 0
print(f"{k} clusters")
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
plt.show()
# Закомментированный код далее - построение отдельных частей графика выше, т.к. исполняется алгоритм ~ 3 минуты
# fig = plt.figure(figsize=(30,10))
# j = 0
# for i in combinations(range(7), 2):
# plt.subplot(3, 7, j + 1)
# plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label2, cmap='viridis', lw=0, s=30)
# plt.xlabel(col[i[1]])
# plt.ylabel(col[i[0]])
# j += 1
# pam = PAM(3)
# pam.fit(data_pca)
# label3 = pam.labels_
# fig = plt.figure(figsize=(30,10))
# j = 0
# for i in combinations(range(7), 2):
# plt.subplot(3, 7, j + 1)
# plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label3, cmap='viridis', lw=0, s=30)
# plt.xlabel(col[i[1]])
# plt.ylabel(col[i[0]])
# plt.legend(label)
# j += 1
# pam = PAM(5)
# pam.fit(data_pca)
# label5 = pam.labels_
# fig = plt.figure(figsize=(30,10))
# j = 0
# for i in combinations(range(7), 2):
# plt.subplot(3, 7, j + 1)
# plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=label5, cmap='viridis', lw=0, s=30)
# plt.xlabel(col[i[1]])
# plt.ylabel(col[i[0]])
# j += 1
import hdbscan
#conda install -c conda-forge hdbscan
hd = hdbscan.HDBSCAN(min_cluster_size=2)
labels = hd.fit_predict(data_pca)
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=labels, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
hd = hdbscan.HDBSCAN(min_cluster_size=3)
labels = hd.fit_predict(data_pca)
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=labels, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
hd = hdbscan.HDBSCAN(min_cluster_size=6)
labels = hd.fit_predict(data_pca)
fig = plt.figure(figsize=(30,10))
j = 0
for i in combinations(range(7), 2):
plt.subplot(3, 7, j + 1)
plt.scatter(data_pca[:, i[0]], data_pca[:, i[1]], c=labels, cmap='viridis', lw=0, s=30)
plt.xlabel(col[i[1]])
plt.ylabel(col[i[0]])
j += 1
km = KMeans(n_clusters=2, random_state=123).fit(data_pca)
label = km.labels_
cl_cent = km.cluster_centers_
ones = np.where(label==1)
zeros = np.where(label==0)
df2_1 = df.loc[df["index"].isin(ones[0])]
df2_0 = df.loc[df["index"].isin(zeros[0])]
df2_0.describe()
df2_1.describe()
km = KMeans(n_clusters=3, random_state=123).fit(data_pca)
label = km.labels_
cl_cent_3 = km.cluster_centers_
twos = np.where(label==2)
ones = np.where(label==1)
zeros = np.where(label==0)
df3_2 = df.loc[df["index"].isin(twos[0])]
df3_1 = df.loc[df["index"].isin(ones[0])]
df3_0 = df.loc[df["index"].isin(zeros[0])]
df3_0.describe()
df3_1.describe()
df3_2.describe()
# conda install -c conda-forge folium
import folium
lat20 = np.array(df2_0["Lat"])
long20 = np.array(df2_0["Long"])
m_cl_0 = folium.Map(location=[df2_0["Lat"].mean(), df2_0["Long"].mean()])
for i in range(len(df2_0)):
folium.Marker([lat20[i], long20[i]], popup=f'The {i} observation').add_to(m_cl_0)
m_cl_0.add_child(folium.ClickForMarker(popup='Waypoint'))
m_cl_0
lat21 = np.array(df2_1["Lat"])
long21 = np.array(df2_1["Long"])
m_cl_1 = folium.Map(location=[df2_1["Lat"].mean(), df2_1["Long"].mean()])
for i in range(len(df2_1)):
folium.Marker([lat21[i], long21[i]], popup=f'The {i} observation').add_to(m_cl_1)
m_cl_1.add_child(folium.ClickForMarker(popup='Waypoint'))
m_cl_1
lat30 = np.array(df3_0["Lat"])
long30 = np.array(df3_0["Long"])
m_cl3_0 = folium.Map(location=[df3_0["Lat"].mean(), df3_0["Long"].mean()])
for i in range(len(df3_0)):
folium.Marker([lat30[i], long30[i]], popup=f'The {i} observation').add_to(m_cl3_0)
m_cl3_0.add_child(folium.ClickForMarker(popup='Waypoint'))
m_cl3_0
lat31 = np.array(df3_1["Lat"])
long31 = np.array(df3_1["Long"])
m_cl3_1 = folium.Map(location=[df3_1["Lat"].mean(), df3_1["Long"].mean()])
for i in range(len(df3_1)):
folium.Marker([lat31[i], long31[i]], popup=f'The {i} observation').add_to(m_cl3_1)
m_cl3_1.add_child(folium.ClickForMarker(popup='Waypoint'))
m_cl3_1
lat32 = np.array(df3_2["Lat"])
long32 = np.array(df3_2["Long"])
m_cl3_2 = folium.Map(location=[df3_2["Lat"].mean(), df3_2["Long"].mean()])
for i in range(len(df3_2)):
folium.Marker([lat32[i], long32[i]], popup=f'The {i} observation').add_to(m_cl3_2)
m_cl3_2.add_child(folium.ClickForMarker(popup='Waypoint'))
m_cl3_2